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ILLINOIS CODE 123~<ftJ. 

Polynomial Approximation 

Complete Program 

About one part in 10 18 

Roughly 20 + [2 .♦ (n/k)]* seconds 

60 

This program prints the coefficients of a polynomial 
which approximates to the function k f (ax ♦ b) with 
error not greater than e over the range c < x < d. 



. The function f is defined by the coefficients 



f p •■ • • ot the Taylor Series . 



£ o> 



t(y) 



+ V .♦ * 2 y* * ,... 



These coefficients may be either given numeri cally or 
generated by an auxiliary subroutine. Any number of 
coefficients up through f 1QO y 100 may be given. 

Code 87 (1.7 Precision Floating Binary Arithmetic) 
is used throughout. 

Machine stops after reading program tape. Then 
insert data tape (see below) and raise black switch. 
There is a sum check at the end of the program tape. 

(1) Enough Taylor coefficients must be given to define 
the function with at least the required accuracy . The. 
residue must remain small compared with e/k for values 
of |y| up to max, ( | ac ♦ b| , |ad ♦ b| ). if in 
doubt, give more coefficients. 

(2) The first number to be printed (to 5 places) is an 
upper limit to the error in the approximation (this 
error may not be reached anywhere in the range, but the 
actual error will not be much smaller). The program 
makes this number as large as- possible without exceeding 
e, by choosing the order of the approximating polynomial . 

This number is followed by the coefficients of the 
polynomial, beginning with the constant term, each to 



DATA TAPE 
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18 places. Numbers are printed in floating decimal 
form with exponent last. 
(3) Remember that the error given is the truncation 

error in replacing an infinite series by a polynomial 

In applying the polynomial, rounding errors will also 

arise . 



! \ 



This program has been found to suffer from a 



serious accumulation of rounding errors if 

(1/20) < (c/d) < 20 
and should not he used in such cases. 
(5) There is no internal check in this program. 
Results should be checked by computing one or two spot 
values of the function. 

Punch first the parameters k, a, b, c, d, e in that 
order, in the form specified for Code Ah (sign, decimal 
digits up to 23; sign of exponent, three 'digit exponent}-. 
Then punch n, the order of the highest Taylor coefficient 
given (n+1 coefficients altogether), with L.H. zeros 
omitted, and followed by the symbol F. Then'- proceed ar. 
follows: 

(a) If Taylor coefficients are given numerically . 
Punch 26 389N followed by the coefficients, punched 
as for Code kk, beginning with the constant term, or 

(b) If Taylor coefficients are generated by auxiliary 
subroutine. 

Punch directive (e.g. 00 IbK) to locate first v rd cr 
auxiliary subroutine, followed by the subroutine, 
followed by a directive (e.g. 26 l6N) to transfer con- 
trol to the subroutine. 
AUXILIARY SUBROUTINE (if used). 

Locations up through 395 are available to hold or be 
used by the auxiliary subroutine. It should place the 
Taylor coefficients as follows: 
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f in 1022,-23 



f in 1020 -21 - in floating binary form, 

f in 1018.-19 etc, I 

and finally transfer control by the order 22 396F 

The following preset parameters may be used, having been set by 
the program tape: 

S3* working positions for Code 87 
Sb* start of Code 87 
S6: location of k 

a, b, c, d, • are stored in 2S6, US6, 6S6, 8S6, 10S6« 

Other preset parameters may be set and used if desired* 

Code 87 may be used. Further numbers may be read (using Code 8?)| 
if so, they should be punched at the end of the tape* 

METHOD 

•1 
The function T (x) - cos (r cos x) can be written as a polynomial 

♦* 
of degree rinx, known as a Tchebyscheff polynomial* The maximum absolute 

value of this function in the range -1 <. x £ +1 is lj of all polynomials of 

degree r satisfying this condition, T (x) has the largest coefficient of x 



r 



In this routine we are concerned with the range c^x^d so we take 

T r (x) - cos[r cos" 1 [(2x - c - d)/(d-c)]] 

which is obtained simply by a linear transformation of the argument • 

Suppose we have a polynomial of degree n (or n+1 terms of a Taylor 

series), 

n 

P n W - ^" *p *?, 

p«0 

and consider the polynomial 

p'W - p n (x)- <vW*>» 

where t^ is the coefficient of x n in T (x)* Clearly P*(x) is of degree n~l§ 

moreover it differs from ?(x) (over the range considered) by at most (a A )<> 

n rr nn' 

If this quantity is less than e, the permitted tolerance^ we can thus reduce 
the degree of the polynomial by 1, Using T n(x) it may then be possible to 
reduce the degree still further* in fact we can proceed until the sum of the 
moduli of the coefficients of the Tchebyscheff s would become greater than e 



The Tchebyscheffs are generated by using the recurrence relation 

T r+1 * "A * Vl 

starting with T Q and T 1 and proceeding to T n . It is only necessary to hold 
two polynomials in the store simultaneously (apart from T,)$ hence two sets of 
locations are reserved to hold coefficients of terms in even and odd Tchebyscheff s 
respectively. 

When the Tchebyscheffs are being used to modify the given polynomial 
it is necessary to re-derive lower order Tchebyscheffs by the relation 

Vi • "A " Vr 

Owing to the way the Tchebyscheffs are. stored, this is done by exactly the same 
orders that generate the Tchebyscheffs • 

The main routine that carries out the above operations is preceded by 
a preliminary routine which makes the shift of origin and change of scale im- 
plied in the parameters k, a and b« 

END COBRECTKW 

The process noimally ends when a further reduction of degree would 
make the sum of the moduli of the coefficients of the Tchebyscheffs greater than 
e$ the current value of that sum is then given as an upper bound to the error* 
This is not quite the most efficient process, for the following reasons . The 
actual error function is the sum of the Tchebyscheffs removed,, The peaks of 
these Tchebyscheffs do not all coincide, so that the peaks of the error function 
will vary in height and will in general be smaller than the quoted maximum error 
(they cannot be. greater). The ideal error function would have all its peaks 
the same height. 

In practice this effect is not often very important because the error 
function is dominated by the last Tchebyscheff to be removed} the others are 
usually much smaller and represent only minor deviations from it. Hence the 
peaks of the error function do not normally differ much in height. Moreover the 
ideal error function can be achieved only by very laborious computation. However^ 
a slight improvement can sometimes be made quite simply, and if this enables a 
further reduction of iegree to be made, this code does it. 

It is assumed that the Tchebyscheff second in importance is that next 
higher in degree to the "dominant" (lowest degree) Tchebyscheff removed. This 
correction can be made only if the coefficient of the second Tchebyscheff is 



-*■ 



is less than half that of the dominant one, in modulus. It aims at reducing 
the effect of the second Tchebyscheff on the peaks of the dominant (similar 
corrections could be made for the other Tchebyscheffs in the error function, 
but these are less likely to be important and are more difficult to program). 

Let c r T p be the dominant Tchebyscheff and c^T x the next 
higher Tchebyscheff in the error function. Simple theor^ yield* a maxUnum 
error of J C r | ♦ | C^J arising frost these, Suppose we now add -c -T 
to the error function. This alters the approximating polynomial but'doer* 
not lower its degree. However, the error function (from these sources) is 
now 

C r T r * Wr + 1 " Wr-1 

c r cos (r cos^x) ♦ c r+1 [cos ((r+Dcos^x) - cos ((r-1 Joes* 3 *)] 

c r cos(r cos" 3 *) - 2Crfl /(l - x*) sin (r cos^x). 

The peaks in the doadnant occur when | cos (r cos" 1 !)] - 1, For such values of 
x, sinfr cos" x) • Oj hence the second term has only a second order effect on 
the height of the peaks. To obtain an upper bound to the error, replace 
/(I - x 2 ) by 1$ this yields the value 

^ C r 2 + ^r*l 2 ) ^ h r i ♦ ^ r ^/\ y . 

The latter is now taken as the upper bound from these sources. 

Before a Tchebyscheff is removed from the given polynomial, a test 
is made to sea whether, on the simple theory, the error would become too large 
If it would not, the Tchebyscheff is removed and the test is repeated on the 
Tchebyscheff of next lower degree. If it would, a further test is made to see 
whether the correction described above would bring the error within the prescribed 
limit. 
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Dec teal Order Input 
00 3K 



8 



00 F 
00 793F 
00 F 
00 £lUF 
00 F 
00 F 
00 F 
00 80gF 
00 F 
00 817F 
00 T 
00 210F 
OOF 
00 UII4F 
00 383K 
00 F 
SO L 
26 SU 
8K 10F 
08 10S6 
8F 2F 
8J UL 
00 F 
IS 6L 
\£ 1016F 
Ul.lF 
26 1009F 
00 39F 
U 2S7 
00 21F 
m 28L 
U6 26L 
SO 8L 



From 13 



n 



work space for Code 8? 



Code 87 



parameters 

constants etc* 
constant term of even 
Tchebyseheffs 
constant term of odd 
Tchebyseheffs. 



Preset 
Parameters 



read parameters 



Use D.O.I, to read n into 2S7 



Clear Q 

Set constant for counting 

in 26L 
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ORDER 




NOTES 


„ 


1 


9 


26 SU 
8L 10L 




waste 






10 


88 (1022)1 


' ?y 12 








11 


8J 11L 

L5 10L 
LO S7 






Read 
coefficients 




12 


kS 10L 
W 26L 
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32 8L 










Ik 


00 39F 
11 J2S7 
00 IF 




Clear Q 






15 


III 28L 
ij2 22L 










16 


OU (20)F 
U6 3S7 


^17« 


Beoa&es Ok (1023 -2n)F 






17 


lh US7 
U6 16L 










18 


L5 2UL 




■■ 






19 


10 S7 
iiO 22tL 
ID 16L 




Change addresses 
and count 






20 


U6 231. 
32'21L 






Shift origin 




21 


22 26L 










■ 


$0 21L 


F*m 20" 








22 


26 SU 










23 


85 ( )F 
8K ( )F 


By 15« 
E&r 20 


Becanes 85(1022-2h)F 






2U 


87 liS6 
OU (102U)F 
OS (102U)F 


fy 19 


inner cycle 






2$ 


8F 2F 
8J 18L 










26 

i 


II 88 ( )F 
50 26L 


Or 8 
From 21 


Becomes 88 (1022 - 2n)F 
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1 
LOCATICW 


OKDEH 


NOTES 


27* 


26 Sb 

8* S6 




k» r 


Scale factors in 


28 


87 (2022)P 
8S (2022)F 


^r 33 




argument and 
function 


29 


& 321 










8$S6 




k» r 




30 


87 2S6 
8S 36 




k.^ x 




31 


8L 281 










00 F 








32 


16 28L 


Tram 29 


Change addresses 






10 S7 




and 


_ 


33 


ijO 28Z, 
10 3S7 


■ 


"" Count 




3U 


36 38SU 
50 2S7 


: 


Set H(liS7> 


- 


3$ 


00 6QP 




■ «— 




k6 US7 




• 8* 2nP 

rtS» 




001*19* 







1 


(22 L) 
(1*9 C8)F 

70 F 
70 F 


-pSF'2 1 

{From 3 


Waste *"" 

Padding 
for EAR 


Clear apace for 
Tchebyscheffs j 


2 


FJ L 
li2 L 






. 


3 


32 L 

#> 3L 




— 




U 


26 SU 








8$ 8S6 




Set T Q and T,« 


5 


80 6S6 








8S F 






6 


86F 




Set H{0) - -<o*d)/2 




8S S8 




*(2) - h/(*t) 


7 


8US8 








8S UF 




1(6) - e (regaining tolerance) 
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MOTES 



8 


86F 




81 b09bS9 


9 


83 b09bS9 




8S 2F 


10 


81 6S6 




80 8S6 


11 


86 2iP 




88 P 


12 


86 2j09l*9 




88 S9 


13 


85 10S6 




8S 6P 


lb 


8K Cb)F 




85 F 


15 


07 S9 




Ob 2S9 


16 


87 2F 




00 S8 


17 


OS S8 




8F 2F 


18 


8J 2bL 




8t 191 


19 


8K ( )F 




85 F 


20 


07 S8 




Ob 2S8 


21 , 


87 2F 




00 S9 


22 


OS S9 




8F 2F 


23 


8J 39L 




8L IhL 


2b 


IS lbL 




LO bS7 


25 


32 38L 



Froa 23 1 

and 69' 
9F bl» 



Font aran Tchabyachaff 



HFram 6l« 
E&r 26* 
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Waste 



- Font odd Tchsfcgraohaff 



V\ S7 



t 
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26 


Ik 1U 


Froa 71 • 


Change degree of odd T and go to form it a 




U6 19L 




- 


27 


U6 bft 








IG 66£ 




■ 


28 


26 38SU 




4 




»> 28L 


Froa 2$, 7k [ 


. 


29 


26 Sk 








_8L (30)L 


Sr 79' 


(lfcate) 


30 


aor 


ftr U2 


Form even coefficient 




01 1022F 




31 


06 S8 








8S US8 




4 


32 


8$ 6F 








82 I4S8 




? 0,K* to use even T 


33 


8S 2S7 








' 83 S7L 






3U 


81 US9 




. 




8S US9 




T 0*1 • if end correction is applied 


35 


86 I4S8 








8S UP 




— 


36 


86 UF 

82 UF 






37 


82 UF 








83 5UL 






38 


8J 82L 




Go to print 




00 F 




Waste 


39 


IS 19L 


From 23 






LO US7 




? use odd T 


IjO 


32 U3L 








16 S7 






Ul 


IU 19L 


From 72 • 


Change degree of even T and go to form it« 




u6 Hii 






1*2 

• 


U6 30L j 
U6 £8L J 






*s ! 


26 38SU ; 


1 






50 U3L j 


From ljO, 76 T 


» 
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Ik- 


26 Sk 








8L (li5)L 


By IV 




w 


8K { )F 
01 1022F* 


9jr 27 




16 


06 S9 




- Form odd coefficient 




8S I4S9 






1*7 


85 6F 








82 I4S9 




- ? O.K« to use odd T 


1*8 


8S 2S7 








83 65L 






Ii9 


81 I4S8 


■v. 




8S I4S8 






5& 


86 I4S9 




? 0*K* if end correction is applied 




8S I4F 






51 


86 UF 
82 I4F 






52 


82 IiF 




- 




83 62L 




~ 


53 


8J 83L 




Go to print 




00 F •" 




Waste 


5U 


87 US9 
8S I4F 


From 37* 




55 


81 2S7 

82 IiF 




- ? correction 0„K. 


56 


83 38L 








8J 77L 




Prepare for end 


57 


8S 6F 


From 33' 






8L 58L 


|~ From 29 » 


Waste 


58 


8K ( )F 
85 I4S8 


-[_ ^r U2« 




59 


07 S8 








OU 1022F 




- Use even T 


60 


OS 1022F 








8F 2F 






61 


8J 70L 







-J 



-7* 



1 

• 

LOCATION 


omm 


wyrtes 


'...,- ;.-■"' 


8t I9L 






62 


87 US8 
85 UP 


Fro* 52* 




63 


81 2S7 




- ? correction 0*1* 


- 


82 liF 






6k 


83 53L 








8j 79L 




Prepare for end 


65 


8S 6F 


From U8« 






8L 66L 


prom Ui' 


Waste 


66 


8K( }F 
85 1^9 


;{j*P % 




67 


07 S9 
om022F 




- Use odd T 


68 


OS 1022F 








8F 2F 






69 


8J 71L 






_ ■■■■'. 


8L 1UL 






70 


LI S7 
22 71L 


From 61 




71 


(26 73L) 


From 69 


Becomes 26 72L, 26 83L $r 73 , ,75S80* 




(26 751) 


From 70' 


Becomes 26261, 26 82L 


72 


■XI S7 


From 71 






26 U1L 






73 


L5 92L 


From 71 


Set link to form lower T ! s 




liO 71L 






7U 


U9 iiS7 




Immobilize T count 




22 28L 






75 


L5 92L 


From 71" 


_ Set link to form lover T's 




10 7U 






76 


1*9 1|S7- 




Imaobilise T count 




22 U3L 






77 


t5 65L 


From 56* 


Set to skip formation of next odd coefficien 




lt2 ua 






78 - 


L5 93L 








22 80L 




' 
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! WCATMN 


ORDER 


Nans 


79 


L5 57L 
li2 29t 


From. 6k* 


Set to skip formation of next even coeff „ 

■ 


80 


L5 9ltL 






81 


ItO 711 
OL 1022 P 
26 38SU 


Proa 78* 


Sit to •ad after next T 


82 


16 12«L 
22 83t 


From 38, 71 


Bad after Forming eiren T 


83 


16 19L 
10 81L 


Prom 53, 71 


lad after forming odd T 


* 


itO L 

56 810, 




Set oount constant for printing 


85 


26 SU 






86 


85 1066 
82 6F 


1 


pTorm maximum error 




89 5F 




Print error 


87 


85 (1022 )F 

89 18f 


By 90 

Prom 88* 


Print terms 


88 


8J 891 






■ 


81 87L 


J- 


89 


16*71 


From 88 






W S7 






90 


U6 87L 




Change address in 87 and count 


91 


36 38SU 




• 


92 


OF P 
26 72L 
26 26L 




Stop 


93 


26 831 
26 26L 






9k 


26 72L 








26 82L 








00 51I4E 








CCDI 87 
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ORDER 
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00 817K 






OS? 


00 2F 






■ 


00 2P 






1S7 


88 F 




' 


! 


OOF 






237 


OOF 








00 { )F 




B9com«3 00 riF 


3S7 


*7< yp 




Bftcoaes 8? (3022 - 2a)F 




00 F 






US7 


m if) 




Becomes 8K 2«F$ hP F 




jl ©OF 




(used to count T's) 




00 16K 




• 




com 108 




Sin check 


■ 

■ 
■ 


20 3838 


1 




. 
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